# CLEAN WORKSPACE AND LOAD PACKAGES --------------------------------------------

rm(list = ls())

library(spmirt)
library(ggplot2)
library(datasim)
library(tidyverse)

# SIMULATE DATA ----------------------------------------------------------------

n <- 300

f <- list(
  prob ~ I(0) +
    gp(list(s1, s2), cor.model = "exp_cor", cor.params = list(phi = 0.04),
       sigma2 = 1),
  size ~ I(1)
  )

data <- sim_model(formula = f,
                  link_inv = list(pnorm, identity),
                  generator = rbinom,
                  n = n,
                  seed = 2
                  )
data <- dplyr::rename(data, gp = gp.list.prob)

ggplot(data, aes(s1, s2)) +
  geom_point(aes(col = gp, size = gp)) +
  scale_colour_distiller(palette = "RdYlBu")

ggplot(data, aes(s1, s2)) +
  geom_point(aes(col = factor(response)), size = 2)

data$gp2 = data$gp + rnorm(300)


vg <- gstat::variogram(gp2 ~ 1, ~ s1 + s2, data, cutoff = 0.7, width = 0.005)
ggplot(vg, aes(dist, gamma, weight = np)) +
  geom_point(aes(size = np)) +
  geom_smooth() +
  expand_limits(y = 0, x = 0) +
  scale_x_continuous(limits = c(0, 0.7))
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

vg <- gstat::variogram(response ~ 1, ~ s1 + s2, data, cutoff = 0.7, width = 0.005)
ggplot(vg, aes(dist, gamma)) +
  geom_point(aes(size = np)) +
  geom_smooth() +
  expand_limits(y = 0, x = 0) +
  scale_x_continuous(limits = c(0, 0.7))
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

noise <- rnorm(n)

# RUN MODELS -------------------------------------------------------------------

Rcpp::sourceCpp("../src/gp-gibss-adap.cpp")
source("../R/ggplot-mcmc.R")

# 50000 -> 17 minutes

# First naive adaptive
iter <- 5 * 10 ^ 4
dist <- as.matrix(dist(dplyr::select(data, s1, s2)))
sigma_prop <- matrix(c(0.138, -0.023, -0.023, 0.1), 2) * 2.38 ^ 2 / 2

system.time(
  samples <- probit_gp_adap(data$response, dist, c(log(1), log(0.04)), iter,
                            sigma_prop)
)
##    user  system elapsed 
## 665.973 509.216 293.875
samples_tib <- as_tibble.spmirt.list(samples, iter/2, select = "param")
summary(samples_tib)
## # A tibble: 2 x 6
##   Parameters `2.5%`  `10%`  `50%`  `90%` `97.5%`
##   <chr>       <dbl>  <dbl>  <dbl>  <dbl>   <dbl>
## 1 V1         0.539  0.687  1.10   1.77    2.31  
## 2 V2         0.0183 0.0226 0.0342 0.0505  0.0623
samples_long <- gather(samples_tib)

as_tibble.spmirt.list(samples, 0, 1, "param") %>%
  gg_trace(alpha = 0.6, wrap = TRUE)

as_tibble.spmirt.list(samples, iter/2, 1, "param") %>%
  gg_density(alpha = 0.5)

as_tibble.spmirt.list(samples, iter/2, 15, "param") %>%
  gg_density2d(V1, V2)

as_tibble.spmirt.list(samples, iter/2, 15, "param") %>%
  gg_density2d(log(V1), log(V2))

nrow(unique(samples_tib)) / nrow(samples_tib)
## [1] 0.1924
coda::effectiveSize(log(samples_tib))
##        V1        V2 
##  433.6714 1068.1642
acf(samples_tib)

cov(log(samples_tib))
##             V1          V2
## V1  0.13624702 -0.01939658
## V2 -0.01939658  0.09665395
# Adaptive with stochastic approximation series

system.time(
  samples2 <- probit_gp_am(data$response, dist, c(log(1), log(0.04)), iter,
                            sigma_prop)
)
##    user  system elapsed 
## 543.614 366.832 227.698
samples_tib2 <- as_tibble.spmirt.list(samples2, iter/2, select = "param")
summary(samples_tib2)
## # A tibble: 2 x 6
##   Parameters `2.5%`  `10%`  `50%`  `90%` `97.5%`
##   <chr>       <dbl>  <dbl>  <dbl>  <dbl>   <dbl>
## 1 V1         0.540  0.699  1.08   1.78    2.31  
## 2 V2         0.0181 0.0226 0.0342 0.0506  0.0634
samples_long2 <- gather(samples_tib2)

as_tibble.spmirt.list(samples2, 0, 1, "param") %>%
  gg_trace(alpha = 0.6, wrap = TRUE)

as_tibble.spmirt.list(samples2, iter/2, 1, "param") %>%
  gg_density(alpha = 0.5)

as_tibble.spmirt.list(samples2, iter/2, 15, "param") %>%
  gg_density2d(V1, V2)

as_tibble.spmirt.list(samples2, iter/2, 15, "param") %>%
  gg_density2d(log(V1), log(V2))

nrow(unique(samples_tib2)) / nrow(samples_tib2)
## [1] 0.16828
coda::effectiveSize(log(samples_tib2))
##        V1        V2 
##  423.8379 1083.1148
acf(samples_tib2)

cov(log(samples_tib2))
##             V1          V2
## V1  0.13389773 -0.02155069
## V2 -0.02155069  0.10361964
# Adaptive with stochastic approximation series

system.time(
  samples3 <- probit_gp_am_scale(data$response, dist, c(log(1), log(0.04)), iter,
                            sigma_prop, 0.35)
)
##    user  system elapsed 
## 548.787 373.370 230.650
samples_tib3 <- as_tibble.spmirt.list(samples3, iter/2, select = "param")
summary(samples_tib3)
## # A tibble: 2 x 6
##   Parameters `2.5%`  `10%`  `50%`  `90%` `97.5%`
##   <chr>       <dbl>  <dbl>  <dbl>  <dbl>   <dbl>
## 1 V1         0.531  0.680  1.08   1.76    2.23  
## 2 V2         0.0182 0.0228 0.0342 0.0509  0.0628
samples_long3 <- gather(samples_tib3)

as_tibble.spmirt.list(samples3, 0, 1, "param") %>%
  gg_trace(alpha = 0.6, wrap = TRUE)

as_tibble.spmirt.list(samples3, iter/2, 1, "param") %>%
  gg_density(alpha = 0.5)

as_tibble.spmirt.list(samples3, iter/2, 15, "param") %>%
  gg_density2d(V1, V2)

as_tibble.spmirt.list(samples3, iter/2, 15, "param") %>%
  gg_density2d(log(V1), log(V2))

nrow(unique(samples_tib3)) / nrow(samples_tib3)
## [1] 0.34392
coda::effectiveSize(log(samples_tib3))
##        V1        V2 
##  684.2511 1249.9670
acf(samples_tib3)

cov(log(samples_tib3))
##            V1          V2
## V1  0.1336862 -0.01520090
## V2 -0.0152009  0.09656268
#        V1        V2
#  671.1306 1195.5363




# system.time(
#   out <- probit_gp_chol2(data$response, dist, c(log(1), log(0.05)), iter, sigma_prop)
# )
# system.time(
#   out0 <- probit_gp(data$response, dist, c(log(1), log(0.05)), iter, sigma_prop)
# )

# plot(out$param[, 2])
# abline(h = 0.03, col = 2)



# # Compare priors with posteriors
# # phi
# hist(exp(rnorm(50000, log(0.03), 0.4)), 200, freq = FALSE )
#

# system.time(
#   out <- probit_gp_chol2(data$response, dist, c(log(1), log(0.05)), iter, sigma_prop)
# )
# system.time(
#   out0 <- probit_gp(data$response, dist, c(log(1), log(0.05)), iter, sigma_prop)
# )

# plot(out$param[, 2])
# abline(h = 0.03, col = 2)



# # Compare priors with posteriors
# # phi
# hist(exp(rnorm(50000, log(0.03), 0.4)), 200, freq = FALSE )
# lines(density(out$param[(iter/2):iter, 2]), col = 2)
# quantile(out$param[(iter/2):iter, 2], c(0.025, 0.975))
#
# # sigma^2
# hist(exp(rnorm(50000, log(1), 0.4)), 200, freq = FALSE, xlim = c(0, 5))
# # hist(exp(rnorm(5000, log(1), 0.3)), 100, freq = FALSE)
# lines(density(out$param[(iter/2):iter, 1]), col = 2)
# quantile(out$param[(iter/2):iter, 1], c(0.025, 0.975))
#
# # phi
# hist(exp(rnorm(50000, log(0.03), 0.4)), 200, freq = FALSE )
# lines(density(out$param[(iter/2):iter, 2]), col = 2)
# quantile(out$param[(iter/2):iter, 2], c(0.025, 0.975))
#
# # sigma^2
# hist(exp(rnorm(50000, log(1), 0.4)), 200, freq = FALSE, xlim = c(0, 5))
# # hist(exp(rnorm(5000, log(1), 0.3)), 100, freq = FALSE)
# lines(density(out$param[(iter/2):iter, 1]), col = 2)
# quantile(out$param[(iter/2):iter, 1], c(0.025, 0.975))
#
# acf(out$param[(iter/2):iter, 1])
# acf(out$param[(iter/2):iter, 2])
#
# acf(out$param[seq((iter/2), iter, 10), 1])
# # acf(out$param[seq((iter/2), iter, 10), 2])
#